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Turbulent  Flow  Validation  in  the  Helios  Strand  Solver 

Andrew  M.  Wissink* 

Aaron  J.  Katz* 

Nicholas  K.  Burgess* 


This  paper  will  present  a  validation  of  the  turbulence  modeling  approaches  used  in  the  Helios  code.  He¬ 
lios  is  a  dual-mesh  overset  unstructured/Cartesian  RANS  solver  developed  by  the  US  Army  for  high-fidelity 
simulation  of  rotorcraft  aeromechanics.  The  paper  will  present  validations  of  standard  3D  flowfields  using 
the  current  near-body  unstructured  code  NSU3D  as  well  as  the  new  strand  mesh  solver  that  is  intended  to  be 
adopted  over  the  next  few  years. 


I.  Introduction 

Rotorcraft  computations  are  challenging  due  to  their  dynamic  and  multidisciplinary  nature.  Rotor  blades  experi¬ 
ence  all  the  complex  aerodynamic  flow  conditions  that  complicate  traditional  fixed-wing  CFD  predictions  -  viscous 
effects,  compressibility,  turbulence,  and  stall  -  except  that  rotors  may  experience  all  these  conditions  within  a  sin¬ 
gle  rotor  revolution.  Moreover,  rotorcraft  flowfields  demand  extremely  accurate  resolution  of  the  wake  vortices  over 
relatively  long  distances  because  of  the  importance  of  blade- vortex  interactions  and  fuselage  effects.  On  the  struc¬ 
tural  dynamics  side,  blades  experience  considerable  aero-elastic  effects  with  highly  unsteady  control  loads  transmitted 
through  a  structurally  complex  hub  connection.  The  strong  coupling  between  the  aerodynamic  and  structural  dynam¬ 
ics  simulations  must  also  be  coupled  with  vehicle  flight  dynamics  and  controls  software  to  achieve  a  trimmed  vehicle 
state.  See  the  reviews  by  Strawn  et  al.1  and  Datta  and  Johnson2  for  further  details  on  the  complexities  of  rotorcraft 
modeling. 

The  Helios  software  was  introduced  as  a  rotary-wing  product  of  the  of  the  CREATE- AV  (air  vehicles)  program,3 
sponsored  by  the  Department  of  Defense  High  Performance  Computing  (HPC)  Modernization  Office,  specifically  to 
provide  a  high-fidelity  analysis  capability  to  the  DoD  for  the  acquisition  of  new  rotary-wing  aircraft.  The  dual-mesh 
paradigm  (Fig.  1)  that  is  the  basis  of  the  CFD  aerodynamics  solution  procedure  consists  of  unstructured  meshes  in  the 
near-body  to  capture  viscous  flow  around  complex  geometry,  and  block  structured  Cartesian  grids  in  the  off-body  to 
resolve  the  wake  through  a  combination  of  high-order  algorithms  and  adaptive  mesh  refinement  (AMR).  An  overset 
procedure  facilitates  data  exchange  between  the  two  mesh  types  as  well  as  enables  relative  motion  between  the  mesh 
systems  -  i.e.  the  near-body  unstructured  rotor  meshes  rotate  and  deform  inside  the  stationary  adaptive  Cartesian 
off-body  grids  system.  Rotor  motion,  deformation,  flight  controls  and  trim  operations  are  provided  by  an  external 
comprehensive  analysis  package.  Coordination  of  the  different  codes  is  managed  through  a  lightweight  and  flexible 
Python-based  infrastructure. 

The  latest  version  of  the  software,  Helios  v3  Rainier,4  was  released  to  government  and  industry  users  in  2012. 
It  is  capable  of  simulating  isolated  or  coupled  rotors  and  fuselages,  multiple  rotors,  and  CFD/CSD  coupled  cases. 
The  both  the  current  near-body  unstructured5  and  off-body  Cartesian6  solvers  have  Detached  Eddy  Simulation  (DES) 
turbulence  modeling  capability.  The  paper  will  present  validations  of  Helios  on  the  turbulent  flow  test  cases  as  defined 
by  the  AIAA  Fluid  Dynamics  Technical  committee. 

A  new  near-body  solver7  is  currently  under  development  that  operates  on  automatically-constructed  Strand8  meshes. 
Validations  with  this  solver  will  also  be  presented. 

II.  DES  Turbulence  Modeling 

Earlier  versions  of  Helios  utilized  the  RANS  solver  in  NSU3D  for  the  near-body  with  the  high-order  inviscid  Euler 
solver  in  SAMARC  in  the  off-body.  This  paradigm  is  considered  sufficient  as  long  as  the  flow  is  attached  but  when  it 

*U.S.  Army  Aerodynamics  Development  Directorate  (AMRDEC),  Moffett  Field  CA 
^Mechanical  and  Aerospace  Engineering,  Utah  State  University,  Logan  UT 
^Science  &  Tech.  Corp,  NASA  Ames  Research  Center,  Moffett  Field  CA 
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Figure  1.  Overset  near/off-body  gridding  paradigm  used  in  Helios.  Unstructured  or  curvilinear  grids  to  capture  geometric  features  and  boundary  layer  near 
body  surface,  adaptive  block-structured  Cartesian  grids  to  capture  far-field  flow  features. 


becomes  separated,  which  occurs  on  the  rotor  in  dynamic  stall  or  on  the  backward  side  of  fuselage  protrusions,  RANS 
tends  to  overpredict  the  lift  and  drag.  In  Rainier  the  solvers  have  been  augmented  in  two  ways.  First,  the  Detached 
Eddy  Simulation  (DES)  has  been  enabled  in  NSU3D  in  order  to  better  model  the  wakes  of  separated  flows.  Second, 
laminar  viscous  terms  have  been  enabled  in  the  off-body.  Further  details  of  these  two  updates  along  with  validation 
results  are  given  in  the  following  sub-sections. 

II.A.  Near-body  DES/DDES  Formulation 

Standard  one  and  two-equation  turbulence  models  used  in  RANS  do  a  steady  average  over  all  turbulent  length  scales 
to  predict  the  turbulent  flow  characteristics.  When  the  turbulence  is  generated  under  time-dependent  conditions,  such 
as  separated  flows,  a  sub-grid  scale  Large  Eddy  Simulation  (LES)  formulation  is  a  better  approximation.  LES  resolves 
the  turbulent  length  scales  to  the  degree  allowed  by  the  computational  mesh,  then  switches  to  a  subgrid  scale  model 
for  the  turbulent  length  scales  that  cannot  be  resolved  by  the  mesh.  Since  LES  is  expensive  if  applied  to  the  very  fine 
viscous  grid  scales  near  the  wall,  a  good  compromise  is  the  hybrid  DES  model,  which  applies  standard  RANS  near  the 
wall  and  switches  to  subgrid- scale  LES  depending  on  grid  density.  That  is,  in  regions  where  the  grid  cannot  support 
the  turbulence  length  scale,  generally  near  wall  boundaries,  are  solved  using  traditional  RANS.  Regions  where  the 
grid  is  length  scale  is  sufficient  to  resolve  the  turbulent  length  scale  are  solved  using  LES. 

Implementation  of  DES  involves  adding  the  LES  model  and  changing  the  length  scale  used  in  the  RANS  Spalart- 
Allmaras  turbulence  model.  The  length  scale  can  be  written  as: 

d  =  d  —  fd  max  (0,  d  -  CDES A)  (1) 

where  d  is  the  distance  from  the  wall,  Cdes  =  0.65,  and  A  a  measure  of  the  local  grid  spacing.9  In  NSU3D,  A  is 
set  to  the  maximum  edge  length  touching  a  given  vertex  on  unstructured  mesh.10  If  fy  is  set  to  1,  then  equation  (1) 
becomes: 


d  =  min  (d,CDESA) 


(2) 


The  length  scale  in  this  original  formulation  is  grid  dependent  causing  incorrect  behavior  to  be  observed  for  some 
cases  with  ambiguous  grid  densities.  The  Delayed-DES  (DDES)  model11  was  introduced  to  be  less  sensitive  to  the 
grid.  DDES  computes  fd  by: 


fd  =  l-tanh([8r</]3) 


(3) 
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where  rj  =  t===  2I2  •  vt  and  v  are  the  kinematic  eddy  viscosity  and  the  molecular  viscosity,  respectively.  Uu  are 

V  UijUijfc  d 

the  velocity  gradients.  K  is  the  Karman  constant.  Either  the  DES  or  DDES  options  can  be  invoked  in  NSU3D. 


z 

z 
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(a)  a  =  45°  (b)  a  =  60° 

Figure  2.  Stalled  NACA  0012  wing  at  M!X!  =  0.15,  computed  density  contours  with  DES  turbulence  model. 


The  DES  capability  is  evaluated  for  calculations  about  a  NACA  0012  airfoil  at  M «>  =  0.15  ,Re  =  105  at  45°  and 
60°  angle  of  attack  (AOA).  The  airfoil  is  fully  stalled  at  these  high  angles  of  attack  and  experiences  time-dependent 
shedding  behavior.  Figure  2(a)  and  (b)  shows  the  density  contour  solution  at  the  y  =  0.8  span  wise  location  for  the  45 
and  60  degrees  AOA  cases.  Figure  3  shows  the  computed  drag  coefficient,  both  the  time  dependent  shedding  behavior 
and  the  average  over  the  simulation.  Figure  4  shows  the  computed  average  Cl  and  Cd  with  DES  and  traditional 
RANS  compared  to  the  experimentally  measured  quantities.12  With  DES  turned  on  the  computed  results  more  closely 
approximate  the  measured  experimental  values. 


NACA001 2  at  a  =  45° 


(a)  a  =  45° 

Figure  3.  Time  dependent  and  average  computed  drag  coefficients,  RANS  vs.  DES. 


NACA001 2  at  cc=  60  0 


(b)  a  =  60° 
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(a)  CL 


(b)  CD 


Figure  4.  Measured  and  computed  average  lift  and  drag  coefficients  using  RANS  and  DES. 


II.B.  DES  Off-Body 

The  conservative  form  of  the  compressible  Reynolds  Averaged  Navier-Stokes  (RANS)  equations  describing  the  con¬ 
servation  of  mass,  momentum  and  total  energy  in  three  dimensions  are  given  as: 

^  +  V  •  (fc(u)  -  Fv(u,  Vu))  =  S(u,  Vu)  (4) 

subject  to  the  appropriate  boundary  and  initial  conditions  within  a  domain  Q.  In  this  work  the  RANS  equations  are 
coupled  to  the  one  equation  turbulence  model  of  Spalart  and  Allmaras  (SA  model)18  with  a  modified  production  term 
due  to  Johnson  and  Allmaras,  given  in  reference.7 


II.B.l.  Spalart  Allmaras  Turbulence  model 


In  this  work  two  forms  of  the  SA  turbulence  model  are  employed:  one  corresponding  to  the  traditional  RANS  version 
(SA-RANS)  and  the  other  corresponding  to  the  DES  version  (SA-DES).  Begining  with  the  RANS  version  of  the 
turbulence  model  which  is  given  by: 


f  +  v,p« 

The  production  term  fP  (v)is  given  as: 
where  S  is  given  according  to 


5  = 


1- 1  [V  •  ( (^  +  pv)  Vv)  +  c&2 p Vv  •  Vv]  -  £>  (u) 

(5) 

P(u)=cblSpv 

(6) 

S  T-  S  S  ^  — cV2S 

sfes+c^S)  - 

f  (cV3-2cV2)5-5 

(7) 

Crot min  (o,  a/ (0*0)  Si jS ji^j 

5  = 


V2/v  2 

k  2d2 


(8) 


where  CD  is  the  vorticity  vector  and  Sij  is  the  strain  rate  tensor  (specified  below).  The  wall  distance  function  is  given 
as: 


d  =  min  (Jw,  0.65AxyZ) 


(9) 
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where  is  the  grid  size  and  dw  is  the  distance  to  the  closest  wall.  In  RANS  mode  d  is  simply  set  to  d  =  dw.  Since 
the  off-body  solver  does  not  directly  contain  any  walls  dw  =  oo.  The  constants  and  functions  fn ,  fV2  are  the  same  as 
those  in  reference.18  The  additional  constant  cV3  is  given  as  cV3  =  .9.  The  destruction  term  is  unmodified  from  the 
original  version  in  reference.18 


dd  (u)  —  cWl  p  fw 


(10) 


The  above  represents  a  conservative  version  of  the  turbulence  model  given  in  reference. 
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II. B. 2.  SA-RANS  System 

The  state  vector  and  flux  vectors  including  those  of  the  SA  model  equation  for  three-dimensional  flow  are  explicitly 
given  as: 


u=  < 


/  \ 

p 

(  \ 
p  u 

/  \ 
pv 

/  \ 
pw 

p  u 

pu2+P 

puv 

puw 

pv 

> ,  F/  =  < 

p  vu 

>,  Ecy=< 

pv2  +P 

>,  Fcz=< 

pvw 

pw 

pwu 

pwv 

pw2  +  P 

E, 

u  (Et  +  P) 

v(Et+P) 

w{Et+P) 

.  P*  , 

,  pvu  J 

pvw 

Fv*  =  < 


Lxy 


UTja  +  VTxy  +  WlXZ  +  Cp  ^ 

5(^  +  PV)I 


Fvz 


\ 

0 

\ 

v 

>* 

F/  =  < 

zyy 

^ yz 

WXyx  +  VT-y-y  +  WTj^  H-  Cp 

( JL  _L  PT  ^ 
yPr  '  Prj  j 

I  ^ 

1  dy 

sfa  +  P'O 

dy 

dy 

> 

1 

r 

0 

1 

(ID 


h  s  = 


u^zx  "f  v^zy  H"  w^zz  A-  Cp  {^pr  +  ^ 

^  +  PV)| 


0 

0 

0 


l  cblSpv+  ^  [c^2pVv- Vv]  —cWlpfw  (|) 


where  p  is  the  fluid  density,  ( u  =  (w,  v,  w))  are  the  Cartesian  velocity  components,  P  is  the  fluid  pressure,  Et  is  the  total 
energy,  cp  is  the  specific  heat  at  constant  pressure,  T  is  the  fluid  temperature,  P,  and  P,  T  are  the  Prandtl  and  turbulent 
Prandtl  numbers  respectively  and  x ^  is  the  total  viscous  stress  tensor  including  the  Boussinesq  approximated  Reynolds 
stresses.  Assuming  a  Newtonian  fluid  and  using  the  Boussinesq  approximation  for  the  Reynolds  stresses,  the  viscous 
stress  tensor  takes  the  form  (with  Xi  =  x, y,  z;  i=  1,2,3): 


j  —  2  (/i  )  Sij 


1  /  dui  +  duj  \  1  duk  g 


2  V  dxj  '  dxi  J  3  dxk  lJ 


(12) 


for  i=  1,2,3  and  j  =  1,2,3 

where  ju  is  the  fluid  viscosity  obtained  via  Sutherland’s  law  and  \ij  is  a  turbulent  eddy  viscosity,  which  is  given  by: 


Pt  = 


fn  ~ 


(13) 


( 

V  t1  ) 


=  7.1 
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It  should  be  understood  that  all  quantities  in  the  above  equations  are  the  Reynolds  Averaged  quantities  (the  usual  () 
notation  is  omitted  for  simplicity).  The  pressure  is  obtained  from  the  ideal  gas  equation  of  state  given  as: 


p  =  (y- 1) 


(u2+v2+w2) 


(14) 


where  y  =  1.4  is  the  ratio  of  specific  heats. 


II. B.  3.  Spatial  Discretization 

The  discretization  of  the  turbulence  model  equation  is  a  critical  components  of  a  CFD  solver.  The  turbulence  model 
has  a  significant  impact  on  the  flow  solution  and  therefore  it  must  be  discretized  with  a  method  that  is  both  sufficiently 
accurate  and  robust.  If  the  discretization  is  too  dissipative  then  insufficient  levels  of  eddy  viscosity  are  produced. ?,? 
However,  recent  work7,  ?’ ?’ ?  has  demonstrated  that  discretizing  the  turbulence  model  to  higher  than  first-order  accuracy 
requires  significant  adjustment  to  both  the  continuous  definition  of  the  turbulence  model  and/or  the  discretization, 
convective  flux  function  and  source  term  treatments  in  order  to  obtain  a  robust  solution.  Furthermore,  reference7  has 
demonstrated  that  the  turbulence  model  discrete  solution  is  very  sensitive  to  the  choice  of  convective  discretization. 
Therefore,  for  explicit  clarity,  the  discretization  of  the  convection  and  diffusion  terms  of  the  turbulence  model  are 
detailed  in  this  section.  The  discretization  used  in  this  work  is  less  than  optimal,  since  the  form  of  the  model  employed 
is  only  advantageous  when  using  an  asymptotically  first-order  accurate  discretization.  However,  the  discretization  is 
robust  and  is  an  example  of  what  can  be  considered  common  practice  for  turbulence  model  discretizations.7,7,7’ 18 


II.B.4.  Convection  Terms 


In  order  to  simplify  the  notation,  the  discretization  is  derived  in  one  dimension,  with  the  understanding  that  in  a  three 
dimensional  setting  the  method  is  applied  to  each  direction  independently.  Since  the  current  version  of  SAMARC 
employes  a  high-order  finite-difference  discretization  for  the  convection  terms  of  the  RANS  equations,  the  turbulence 
model  will  also  be  discretized  using  a  finite-difference  discretization.  In  one  spatial  dimension  the  convection  term  of 
the  SA  turbulence  model  equation  is  given  as: 


dpvt) 


(15) 


which  is  discretized  using  an  asymptotically  first-order  upwind  finite-difference  method.  Consider  that  the  domain 
has  been  gridded  and  the  flow  solution  is  defined  as  a  grid  function  such  that  u  «  {ut}.  The  the  first-order  upwind 
finite-difference  of  the  flux  derivative  at  the  grid  node  i  is  given  as: 


3£pv  _  (PV«),-+i/2  ~  (pv»)i-i/2 

dx  Ax 


(16) 


This  is  essentially  a  decoupled  flux  formulation  where  the  turbulence  model  is  upwinded  based  solely  on  the  direction 
of  the  convective  velocity  and  up  winding  due  to  acoustic  propagation  are  ignored.7  In  order  to  maintain  a  first-order 
accurate  discretization  and  also  maintain  pure  upwind  approach  for  a  constant  velocity  field  the  flux  at  the  half  node 
i- b  1/2  is  given  as 

(pvw)i'_j_i/2  =  j  ((pv),  +  (pv)i+1  +sign(ui)  ((pv),.-  (pv)i+1))  (17) 

and  the  flux  at  the  half  node  i  —  1/2  is  given  as: 

(pvw)i_i/2  =  j  ((pv)*_i  +  (P *)i  +  sign(ui)  ((pv),.!  -  (pv),.))  (18) 

In  essence  this  method  convects  the  turbulence  model  variable  pv  into  the  control  volume  surrounding  the  node  i  in  an 
upwind  fashion  at  velocity  Ui.  Note  that  this  approach  defines  the  velocity  at  the  half  nodes  using  a  piecewise  constant 
extrapolation  form  node  i. 


II. B.  5.  Diffusion  Terms 

One  of  the  principal  concerns  with  discretizing  the  SA  turbulence  model  equation  is  the  diffusion  term  treatment.  The 
diffusion  term  consists  of  two  parts  given  as: 


1 

a 


V  •  ((A/  +  pv)Vv)  +  Q,2pVv-Vv 

s - V - -  s - V - ' 

1  2 


(19) 
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While  the  term  labeled  1  in  Eq.  (19)  can  discretized  without  difficultly  there  is  some  debate  as  to  how  to  discretize 
the  second  term.  If  one  were  to  be  general  and  rigorous  one  treat  this  as  a  source  term  since  it  cannot  be  case  in 
a  divergence  form.  However,  since  a  first-order  upwind  formulation  is  employed  for  the  convective  fluxes  one  can 
re-write  the  entire  diffusion  term  as: 

~  [V-  ((p  +  pv)  Vv)  +Q,2pVv  •  Vv]  =  1  [V  •  ((^+  pv  (1  +q,2))  Vv)  —  Q,2pvV2v]  (20) 


which  is  the  so-called  non-conservative  diffusion  term.18  Since  the  solution  is  asymptotically  first-order  accurate  one 
can  regard  the  pv  in  the  second  term  as  constant  on  the  stencil  and  discretize  the  laplacian  of  V  using  a  second-order 
accurate  central  difference  as: 


cb2  pvV2v 


Cb2  (pv),- 


( V,+l  —  2v,  +v,_i  \ 

V  Ajc2  ) 


(21) 


The  first  diffusion  term  is  discretized  using  a  generalization  of  the  central  difference  of  the  second  derivative: 


V.((A/  +  pv(l+C,2))Vv)  W  -E  (- ^'+1  +  \  (v«+l  -2V,-+V,-_i)  +  (V«-l  -V,))  (22) 

m  =  Bi  +  pv*  ( 1  +  Cbl ) 


Adding  these  two  pieces  together  and  grouping  the  terms  by  the  stencil  point  gives  the  following  final  formula  for  the 
discretization  of  the  entire  diffusion  term  as: 


1  [V-((p  +  pv(l+cfo2))Vv)-c&2pvV2v]  W  Al(vi+1D,-+1-vA-(-v,_iA-_i) 


A  =  +  rh  +  - V, 

rp+i+ri/ 

A+i  =  — ^ - cb2pVi 


(23) 


This  formula  is  formally  second-order  accurate,  but  since  the  convection  terms  are  first-order  accurate  the  second-order 
accuracy  is  not  critical. 


II.B.6.  Source  Terms 

The  source  term  discretization  is  trivial  compared  to  the  convection  and  diffusion  terms.  The  source  terms  are  simply 
evaluated  point-wise  for  each  node  in  the  grid.  However,  the  vorticity  and  strain  terms  required  for  the  production 
term  require  gradients  of  the  velocity  field.  These  gradients  are  obtained  via  a  central  difference  of  the  first  derivative: 

du 
dx 

from  which  the  vorticity  and  strain  rate  magnitude  are  derived  using  the  standard  formulas: 


Mi-\- 1  Ui—  l 

2Ax 


(24) 


/->  ^  /  fdw  dv\ 2  ( du  dw\ 2  7 dv 

+U-^j  +(s 

SijSji  =  \Js j ,  +  +  2  o  j  —  Si,  —  2Si.  —  Sfj 

o  _  1  (dm  duj\ 

*iJ  _  2  \dxj  +  dxi  ) 


where  all  the  derivatives  are  evaluated  at  a  grid  node  according  to  the  formula  in  Eq.  (24). 


(25) 


II.B.7.  Explicit  Time-stepping 

The  S  AMARC  solver  was  originally  developed  to  solve  the  Euler  equations  in  the  off-body  region  of  the  computational 
domain,  and  consequently  employs  explicit  time-stepping  for  temporal  discretization.  The  explicit  nature  of  S AMARC 
has  become  deeply  embedding  within  the  design  of  the  software  and  it  is  beyond  the  scope  of  this  work  to  change  this 
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fundamental  assumption.  Therefore  in  order  to  ensure  the  robustness  of  a  turbulent  version  of  SAMARC,  the  time-step 
must  be  redefined  to  account  for  the  additional  complexity  of  the  turbulence  model  source  and  diffusion  terms. 

In  order  to  develop  the  time-step  definition  for  the  turbulent  version  of  SAMARC  a  linear  convection  diffusion 
model  equation  is  employed  to  establish  the  stable  explicit  time-step.  Furthermore,  the  equally  spaced  Cartesian 
nature  of  the  SAMARC  grids  allows  one  to  do  analysis  in  one  spatial  dimension.  Multi-dimensional  results  come 
summing  the  results  in  three  dimensions.  The  linear  convection  diffusion  model  equation  is  give  as: 


du 

dt 


d(au)  d2u 

— — -  —  v - 

dx  dx2 


=  Su 


(26) 


Eq.  (26)  is  discretized  using  a  first  order  upwind  scheme  for  the  convection  term  and  a  second  order  central 
difference  for  the  diffusion  term.  The  time  derivative  is  discretized  using  a  first  order  forward  Euler  method.  The  fully 
discretized  Eq.  (26)  is 


u’l+l  -u’l 
At 


Ax 


7  i  wf,  i  —2 u!J  +  i 

_v_Lhl -  , - ^  =  Sun 


Ac 2 


(27) 


the  equation  for  new  time  level  value  of  n”+1  is 


M”+1  =  m"  —  At  a  - 


Ax 


—  +Atv’^±l- 


-2m? 


h- 1 


Ax2 


+  A  tSifi 


(28) 


The  stability  of  of  this  scheme  is  analyzed  by  employing  VonNeumann  analysis  on  Eq.  (28).  VonNeumann  analysis 
is  valid  for  period  boundary  conditions  and  an  equally  spaced  mesh  containing  N  points  covering  a  domain  of  length 
L.  One  starts  by  introducing  the  finite  Fourier  Transform  of  the  discrete  solution 


N 

2 


*?  =  I 

k=~  f 

i  =  V^T 


'clTikx; 


(29) 


Assuming  a  periodic  solution  and  boundary  conditions  one  can  write  the  coordinate  xj  =  i  •  Ax.  Furthermore,  to  keep 
the  book  keeping  simply  one  can  define  27t^Ax  =  ^.  It  is  sufficient  to  consider  an  arbitrary  Fourier  mode  for  the  stability 
analysis.  Substituting  a  single  Fourier  mode  for  the  solution  into  Eq.  (28)  yields 


—  At  a 


une^  -  une^-^ 


Ax, 


+  Arv 


%neK(i+\)  _  2 une^  +  une V'-1) 
Ax2 


H-  At  Si/1  c 


^(9 


(30) 


Simplification  and  division  by  un  gives 


8  = 


un+l 


Ata 

Ax 


Av 

Ax2 


_2  +  ^(-D) 


T-  At  S 


Defining  the  Cell  Reynolds  number  (Re a*)  and  CFL  number  (a)  as 


(31) 


Re  Ax 


a 


aAx 

v 

Ata 

Ax 


(32) 


g  =  1  -  a  (I  -  e^(_1))  -  -  2  +  +  AS  (33) 

Figure  5(a)  shows  an  example  Fourier  footprint  using  CFL  =  .75  and  Re  a*  =  10  with  zero  source  term,  one 
should  note  that  all  modes  are  stable.  Figure  5(b)  shows  the  Fourier  symbol  footprint  with  a  non  zero  source  term, 
demonstrating  the  onset  of  instability  as  a  result  of  adding  the  source  to  the  equation.  The  magnitude  of  the  Fourier 
symbol  is 


\g\2  =  [(1  —a  —  2a  +  SAt)  +  (a  +  2a)  cos  (^)]2  +  c2  sin2  (£) 

a 

a  =  — — 

Rsax 


(34) 


8  of  25 


American  Institute  of  Aeronautics  and  Astronautics 


11 


Real(g)  Rea  1(g) 

(a)  Fourier  symbol  footprint,  S  =  0  (b)  Fourier  symbol  footprint,  S  =  —  1 

Figure  5.  Example  Fourier  symbol  footprint  plots  using  a  CFL  =  .75  and  Re  a*  =  10,  a  zero  and  non  zero  source  term. 


which  has  extrema  at  £  =  ±71.  Evaluting  Eq.  (34)  at  the  £  =  ±7T  gives 

max(\g\2)  =  (—  1  ±4oc  —  SAf  ±2o)2 


(35) 


The  criteria  for  stability  is  therefore 


i±  <  i 

—  1  ±  4a  —  SAt  ±  2a  ±  ±  1 


Af  < 


4_y_ . 
4  A*2 


o  _ o 

'ZAx  ^ 


(36) 


This  analysis  holds  for  linear  problems  and  before  applying  it  to  determine  the  time-step  for  SAMARC  one  must 
properly  interpret  and  apply  the  results,  in  light  of  the  non-linearity  of  the  Navier-Stokes  equation.  Firstly,  the  value  of 
a  is  taken  as  aoo ,  which  is  the  freesream  speed  of  sound.  Coincidentally  this  is  the  estimate  of  the  maximum  eigenvalue 
currently  used  in  the  inviscid  version  of  SAMRC.  Secondly,  v  which  is  the  viscosity  in  the  analysis  should  represent 
the  maximum  possible  diffusion  coefficient  for  the  SA  turbulence  model. 


v  =iu±max(pv) 

X 


(37) 


Finally,  the  SA  turbulence  model  has  non-linear  source  terms  that  my  take  on  positive  or  negative  values.  The  analysis 
assumed  that  S  was  a  positive  number  therefore  the  value  of  S  for  the  Navier-Stokes  equations  with  SA  turbulence 
model  is  defined  as  the  minmum  of  zero  and  the  linearization  of  the  SA  turbulence  model  source  term  with  respect  to 


pv. 


S  =  min 


0,min 

X 


9(fP(u)  —  2>(u)) 

3u 


(38) 


which  will  ensure  that  source  term  contributions  only  make  the  time- step  smaller.  Putting  all  this  together  and  account¬ 
ing  for  the  fact  that  there  are  3  spatial  derivative  terms  gives  the  final  time-step  restriction  that  is  used  in  SAMARC  for 
turbulent  computations 


At  < 


/tA'+max±pv) 

Ax2 


±2^ 


—  min  ^0,  i 


Y«!h2W) 


5u 


)) 


(39) 


III.  Strand  Mesh  Solver 

The  overset  dual-mesh  “strand’-Cartesian  approach  has  been  proposed  and  studied  in  earlier  works7, 8,13-17  as  a 
viable  means  to  support  automatic  viscous  mesh  generation  and  adaptation.  In  the  strand  paradigm,  a  body-fitted  near¬ 
body  mesh  is  constructed  by  a  set  of  straight  line  segments  grown  directly  from  the  surface,  each  with  the  same  point 
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Figure  6.  Strand- Cartesian  grid  system. 


distribution  in  the  normal  direction,  forming  a  thin  layer  of  mostly  prismatic  elements  around  the  body.  Once  outside 
the  viscous  boundary  layer,  strands  transition  to  isotropic  block  structured  Cartesian  grids.  The  two  grid  systems 
intersect  through  overlapping  chimera  overset  procedures.  The  procedure  is  similar  in  concept  to  standard  prismatic 
unstructured  grid  generation  techniques,  in  which  prismatic  cells  are  grown  at  the  surface  in  the  viscous  boundary 
layer  with  tetrahedra  elsewhere,  except  in  the  strand  approach  Cartesian  grids  are  used  in  place  of  tetrahedra  for  the 
Euler  solution. 

In  addition  to  streamlined  and  automatic  meshing  capability,  the  strand-Cartesian  approach  presents  three  other 
important  advantages.  First,  both  strand  and  Cartesian  meshes  may  be  represented  with  extremely  low  memory  de¬ 
scriptions,  enabling  the  entire  global  mesh  description  to  fit  on  each  processor  in  a  parallel  environment.  This  allows 
for  significant  gains  in  efficiency  and  scalability  of  domain  connectivity,  effectively  eliminating  inter-processor  search 
routines.  The  savings  become  even  more  significant  in  the  case  of  moving  body  simulations  for  which  domain  con¬ 
nectivity  must  be  re-established  at  each  unsteady  time-step.  Second,  both  strand  and  Cartesian  meshes  possess  at  least 
some  grid  structure,  facilitating  efficient  implementations  of  high-order  accurate  discretizations  and  solution  meth¬ 
ods.  These  methods  include  high-order  finite  differencing,  line-implicit  solvers,  and  directional  multigrid  coarsening. 
Third,  both  the  strand  and  Cartesian  grids  easily  permit  use  of  Adaptive  Mesh  Refinement  (AMR).  Because  all  strands 
use  the  same  normal  point  distribution,  adaptation  is  entirely  surface-based.  This  avoids  cell  quality  and  edge  swap¬ 
ping  complexities  that  have  traditionally  plagued  volume-based  unstructured  AMR.  AMR  on  Cartesian  grids  has  been 
known  for  years  to  be  very  effective  because  the  logical  data  structure  naturally  facilitates  a  hierarchical  mesh  repre¬ 
sentation  and  Cartesian  cells  do  not  suffer  cell  quality  issues  with  frequent  and  persistent  adaptation,  as  can  occur  with 
tetrahedral  elements. 

III.A.  Strand  Mesh  Generation 

The  starting  point  for  the  mesh  generation  is  a  tessellated  surface  composed  of  either  triangles,  quadrilaterals,  or  a 
mix  of  both.  Strands  consist  of  straight  line  segments  of  equal  length,  number  of  points,  and  point  distribution,  grown 
from  surface  vertices.  The  strands  initially  are  grown  normal  to  the  surface  and  then  smoothed  to  provide  coverage  in 
convex  corners  (Fig.  7a  &  b)  and  to  push  crossing  strands  away  from  the  surface  in  concave  comers  (Fig.  7c  &  d).  The 
desired  degree  of  smoothing  in  the  mesh  may  be  adjusted  at  runtime  through  an  input  parameter. 

Spacing  along  each  strand  ranges  from  viscous  at  the  root  to  transitional,  or  “Euler”  spacing  8 e  at  the  tip.  The  set 
of  strands  produce  a  prism  stack  associated  with  each  surface  triangle.  Any  negative  volumes  associated  with  crossing 
strands  are  clipped  by  associating  an  integer  “clip  index”  icup  with  each  surface  triangle.  The  clip  index  may  also  be 
used  to  clip  elements  that  protrude  an  outer  mold  line  for  two  strand  meshes  that  lie  in  close  proximity  to  one  another. 
The  user  can  supply  the  desired  strand  length  Ls  through  input  or  allow  the  strand  length  to  be  computed  automatically. 
If  the  length  Ls  is  computed  automatically,  the  algorithm  seeks  to  make  the  transitional  spacing  8 e  at  the  strand  end 
equal  to  the  surface  tessellation  spacing.  This  ensures  the  transition  cells  are  roughly  isotropic  at  the  strand  ends  so 
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' 

(b)  Convex  corner  -  Smoothed 


(c)  Concave  corner  -  Non-smoothed 


(d)  Concave  corner  -  Smoothed 


Figure  7.  Strand  direction  vector  smoothing. 
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they  transition  nicely  to  Cartesian  off-body  meshes.  The  strand  length  Ls ,  Euler  spacing  8 £,  and  clip  index  icup,  are 
all  pictured  graphically  in  Figure  6. 

Once  a  near-body  strand  mesh  is  available,  an  adaptive  Cartesian  off-body  mesh  is  automatically  constructed. 
The  output  of  the  near-body  strand  mesh  generation  is  a  set  of  prism  stacks  each  with  a  designated  clip  index  iciip 
and  the  normal  spacing  of  the  prism  element  at  this  clip  index,  or  Euler  spacing  8 e,  provides  the  basis  for  the  initial 
Cartesian  grid  generation.  The  location  of  the  clip  index  elements  (x,y,z)  and  the  Euler  spacing  8 e  are  provided  to  the 
Cartesian  grid  generator.  Block  structured  Cartesian  grids  are  built  in  a  hierarchical  fashion,  the  coarsest  level  defines 
the  physical  extent  of  the  computational  domain  and  new  levels  are  constructed  from  coarsest  to  finest.  Each  finer  level 
is  formed  by  selecting  cells  on  the  coarser  level  and  then  clustering  the  marked  cells  together  to  form  block  regions 
that  will  constitute  the  new  finer  level.  The  Cartesian  grid  is  initially  refined  to  match  the  Euler  spacing  8 e  elements 
in  the  strand  mesh,  and  Cartesian  grids  are  then  adapted  throughout  the  simulation  to  capture  time-dependent  solution 
features  such  as  vorticity.  Further  details  of  the  off-body  mesh  generation  procedure  have  been  presented  in  previous 
work.17 


IV.  RANS-SA  Solver 


In  this  work  we  solve  the  Reynolds-averaged  Navier-Stokes  (RANS)  equations  in  three  dimensions.  Turbulence 
closure  is  accomplished  with  the  Spalart-Allmaras  (SA)  model.18  The  RANS-SA  equations  may  be  expressed  as 


d_Qd_F1J_Fl  = 

dt  dxj  dxj  ’ 


(40) 


where  the  conserved  variables,  Q ,  inviscid  fluxes,  Fj,  viscous  fluxes,  FJ,  and  source  term,  S ,  are  defined  as 


p ) 

(  P  Uj  \ 
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5  Fj  = 

pUiUj  +  p8/y 
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0 
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p  huj 

5  j 

Gijui  ~  Qj 

0 

w 

V  pvw j  ) 

\  / 

\  O  dxj  / 

(fp-©+cMPgg/ 

Here,  p  is  the  density,  Uj  is  the  Cartesian  velocity  vector,  e  is  the  total  energy  per  unit  mass,  V  is  the  turbulence  working 
variable,  p  is  the  pressure,  h  is  the  total  enthalpy  per  unit  mass,  O/j  is  the  deviatoric  stress  tensor,  qj  is  the  heat  flux 
vector,  and  r|/o  is  the  turbulent  diffusion  coefficient.  The  turbulent  source  term  consists  of  a  production  term,  fP,  and 
a  destruction  term  (D .  The  stress  tensor  is  defined  as 


°ij  =  2(M+IJT)Sij, 


(42) 


where  p  is  the  dynamic  viscosity,  pj  is  the  turbulent  viscosity,  and  Sij  is  the  rate  of  strain  tensor,  defined  as 


1  (  dui  duj  \  1  duk 

Slj  2  V  dx;  )  3  dxi,  u 


The  heat  flux  vector  is  obtained  with  Fourier’s  Law, 


Qj  ~  Cp 


pr 


Pr  Prj  )  dxj  ’ 


dT 


(43) 


(44) 


where  Cp  is  the  specific  heat,  Pr  is  the  Prandtl  number,  Prj  is  the  turbulent  Prandtl  number,  and  T  is  the  temperature. 
The  ideal  gas  equation  of  state,  p  =  p RT  is  used  to  close  the  equations. 


IV.A.  Discretization  and  Solution  Methods 

The  strand  grid  spatial  discretization  is  based  on  a  cell-centered  approach  where  the  primary  unknowns  are  located 
at  the  centroid  of  the  prisms  formed  by  adjacent  strands.  The  solver  accommodates  both  quadrilateral  and  triangular 
prisms  depending  on  the  surface  topology.  However,  control  volumes  are  composed  entirely  of  triangular  facets  by 
triangulating  any  non-planar  quadrilateral  faces.  This  is  important  for  second-order  accuracy  on  general  prismatic  grids 
with  no  assumption  of  underlying  smoothness.19  Linear  reconstruction  is  employed  to  obtain  second-order  accuracy 
through  first  obtaining  consistent  nodal  values  of  the  conserved  variables  from  surrounding  cell-center  values.  A 
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projection  method  is  used  to  obtain  these  nodal  values  via  least  squares  interpolation  in  a  regression  plane  through  the 
three-dimensional  stencil  of  cells  surrounding  a  strand.19  Once  the  nodal  values  have  been  obtained,  a  Green-Gauss 
surface  integration  procedure  is  performed  to  obtain  cell  gradients  in  each  control  volume. 

Inviscid  fluxes  rely  on  a  reconstruction  upwind  formula  for  the  numerical  flux  based  on  the  approximate  Riemann 
solver  of  Roe,20 

f  =  \  (!F(Qr)  +  HQl))  -  l  \A{Qr,Ql)\  (Or  -  Ql ) ,  (45) 

where  J  =  Fjnj  is  the  directed  flux  at  a  face  with  normal  nj,  and  A  =  /dQ  is  the  directed  flux  Jacobian.  The 
viscous  terms  are  computed  using  values  of  Q  and  Vg  determined  at  each  face, 


F  =  r(Qf,VQf), 


(46) 


where  /  refers  to  the  face  reconstructed  values.  These  face  values  are  easily  obtained  once  nodal  values  have  been 
reconstructed  using  the  projection  method  described  above.  This  method  is  similar  to  the  node  averaging  schemes  out¬ 
lined  by  Diskin,  et  al.21  Both  the  inviscid  and  viscous  discretization  methods  described  herein  have  been  verified  to  be 
second-order  accurate  for  arbitrary  prismatic  meshes  under  a  variety  of  conditions  using  the  method  of  manufactured 
solutions.19 

The  result  of  the  spatial  discretization  of  the  viscous  and  inviscid  fluxes  is  a  coupled  set  of  non-linear  equations. 
In  this  work,  we  adopt  a  pseudo-time  framework  to  march  the  steady  or  unsteady  discretized  equations  to  steady-state, 

dQ  ,  x 

V^+R(2)=0.  (47) 


Here,  V  is  the  cell  volume,  and  x^  is  the  pseudo-time  variable.  The  residual,  R(Q),  contains  the  inviscid  and  viscous 
flux  balances  at  each  cell  based  on  the  cell-center  discretization  schemes  described  above.  In  order  to  reach  a  pseudo¬ 
steady  state  using  an  implicit  scheme,  the  residual  is  linearized,  leading  to  the  following  linear  system  to  be  solved  at 
each  pseudo-time  step: 


V  dRk 
Ax/  +  dQ 


- R{Qk )• 


(48) 


Here,  dRk/dQ  is  the  Jacobian  of  the  residual.  The  linear  system  in  Equation  47  in  general  is  large  and  sparse,  rendering 
direct  inversion  impractical.  Iterative  line  Gauss-Seidel  (GS)  methods  are  employed  to  solve  this  system,  where 
contributions  along  strands  are  collected  to  form  a  tridiagonal  system.  To  facilitate  the  line  GS  iterations  and  to 
increase  robustness,  we  introduce  an  additional  “linear  time”  variable,  X/, 


dQ 

VJT 

dx, 


V  dRk 
A  xkI+  dQ 


(Qk+l-Qk)  =  —R(Qk). 


(49) 


The  linear  time  is  introduced  to  improve  the  diagonal  dominance  of  the  line  GS  procedure  in  order  to  increase  robust¬ 
ness.  Rearranging  Equation  49  in  terms  of  solution  updates  in  linear  time  results  in 


VIA- 


dRk 

dQ 


- mk ) 


V  dRk 

- /+  -^r— 

Ax*  dQ 


(Ql-Qk)- 


(50) 


Upon  convergence  of  the  linear  iterations  in  /,  the  linear  system  of  Equation  48  is  satisfied.  At  that  point,  the  next 
pseudo-time  step  in  k  proceeds.  When  the  pseudo-time  iterations  converge,  then  the  residual  equation  R{Q)  is  satisfied 
for  a  given  physical  time  station.  All  Jacobian  terms  in  this  work  are  first  order  and  retain  only  nearest  neighbor 
contributions.  Further  details  of  the  implicit  solution  method  may  be  found  in  previous  work.16 


IV.B.  Turbulence  Model 


The  standard  SA  model  is  used  when  the  turbulent  working  variable  is  positive.  Details  of  the  positive  model,  including 
the  well-known  definitions  of  the  production  and  destruction  terms,  may  be  found  in  the  original  work  by  Spalart  and 
Allmaras.18  Modifications  to  the  model  to  accommodate  negative  values  of  the  turbulence  working  variable  have  been 
suggested  recently  by  Allmaras22  and  are  employed  in  this  work.  In  the  case  of  negative  values  of  V,  the  following 
turbulence  equation  replaces  the  standard  model: 


3v  dv  ( v 

¥+“J'a^  =  CM(1_Q3)n^+CM'1  \d 


d  ( .  _  dv  \  dv  dv 

d^l(V  +  Vfn)d^J+Cb2d^d^ 


(51) 
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where, 

f  Cm+X3  (n  ^ 

Jn=  r  _  3  ,  (C„i  =  16). 

Here,  £2  is  the  vorticity  magnitude,  d  is  the  distance  to  the  nearest  wall,  and  %  =  v/v  is  the  ratio  of  the  turbulent 
working  variable  to  the  kinematic  viscosity  of  the  fluid.  All  other  constants  in  Equation  5 1  take  the  values  found  in 
the  standard  model. 


V.  Results 


V.A.  Helios  Validation 

V.A.  1.  Flow  over  a  sphere 

The  first  application  example  is  the  flow  over  a  sphere  at  high  reynolds  number.  This  is  a  canonical  fluid  mechanics 
problem  and  provides  a  suitable  initial  test  of  the  turbulence  model  in  SAM  ARC.  The  flow  conditions  for  this  case 
are  free-stream  Mach  number  =  .3,  angle  of  attack  a  =  0.0,  side  slip  angle  P  =  0.0  and  a  Reynolds  number  based 
on  sphere  diameter  of  Re  =  6.7 6e6.  The  near-body  mesh  contains  1,440,699  grid  points  and  the  off-body  mesh  has 
139, 485,538  grid  points.  A  slice  of  the  mesh  at  the  y  =  0  plane  is  shown  in  Figure  8(a).  This  test  case  is  time-accurate 


(b)  Close-up  of  near-body  mesh  ,  y  =  0  plane 


Figure  8.  Mesh  used  for  the  computation  of  the  flow  over  a  sphere. 


with  a  time-step  of  1.0e~7  seconds,  a  total  of  50,000  time-steps  are  simulated. 

This  work  is  focused  on  evaluating  the  various  physical  modeling  options  in  SAM  ARC.  Therefore  the  physical 
model  used  in  the  near-body  mesh  is  always  viscous  and  employs  either  the  SA-RANS  or  SA-DDES  turbulence 
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models.  The  overset  dual-mesh  capability  of  Helios  allows  for  various  combinations  of  physical  models  in  the  near¬ 
body  (NB)  and  off-body  (OB)  meshes.  The  following  combinations  of  physical  models  are  employed  for  this  test 
case:  SA-RANS  in  the  near-body  with  inviscid  off-body,  SA-RANS  near-body  with  SA-RANS  off-body,  and  SA- 
DDES  near-body  with  SA-DES  off-body.  Figure  9(a)  depicts  the  vorticity  contours  at  the  final  time-step  on  the 

Vorticity  Magnitude  [PLOT3D] 
o 

0.800 

■  ■  ■ 

0.005 


(a)  Vorticity  Contours  ,y  =  0  plane 

Figure  9.  Vorticity  for  the  flow  over  a  sphere  with  an  inviscid  solver  in  the  off-body  region. 


Vorticity  Magnitude  [PLOT3D] 
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(a)  Vorticity  Magnitude,  y  =  0 

Fi 

o 

2115.563 

E  ■  ■ 

0.000 


(b)  Eddy  Viscosity,  y  =  0  plane 

Figure  10.  Vorticity  for  the  flow  over  a  sphere  with  a  SA-RANS  solver  in  the  off-body  region. 
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y  =  0  plane  for  using  the  inviscid  off-body  and  RANS  near-body  physical  model  combination.  The  inviscid  version  of 
the  off-body  solver  generates  a  highly  unsteady  wake  solution,  which  represents  the  base-line  Helios  solution  for  this 
problem. 

Figure  10(a)  shows  the  vorticity  contours  on  the  y  =  0  plane  using  SA-RANS  near-body  and  SA-RANS  off-body 
combination.  Employing  the  SA-RANS  turbulence  model  in  the  off-body  region  has  changed  the  wake  structure  from 
a  length  scale  rich  and  obviously  unsteady  wake  to  one  that  more  closely  resembles  that  of  a  steady  wake  where  all 
scales  have  been  modeled.  Figure  10(b)  shows  the  turbulent  eddy  viscosity  for  this  case,  which  is  highest  near  the  body 
and  quickly  decreases  in  the  downstream  region.  The  eddy  viscosity  distribution  confirms  that  extensive  turbulence 
modeling  is  occurring  in  this  region. 

Vorticity  Magnitude  [PLOT3D] 
o 

0.800 

■  ■  ■ 

0.005 


(a)  Vorticity  Magnitude,  y  =  0 


FI 
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200.000 


0.000 


(b)  Eddy  Viscosity,  y  =  0  plane 

Figure  11.  Vorticity  for  the  flow  over  a  sphere  with  a  SA-DES  turbulence  model  in  the  off-body  region  and  a  DDES  turbulence  model  in  the  near-body  region. 


Figure  1 1(a)  depicts  the  vorticity  contours  on  the  y  =  0  plane  using  SA-DDES  in  the  near-body  mesh  and  SA-DES 
in  the  off-body  mesh.  Employing  SA-DES,  which  is  partially  a  large  eddy  simulation  (LES)  method,  in  the  off-body 
region  results  in  a  wake  which  is  again  unsteady  and  contains  a  wide  variety  of  captured  length  scales.  Examination 
of  Figure  11(b),  which  depicts  the  eddy  viscosity  using  SA-DES  in  the  off-body,  confirms  that  the  turbulence  model 
is  active  in  the  wake  region.  Comparing  Figure  9(a)  to  Figure  11(a)  qualitatively  demonstrates  that  using  an  inviscid 
solver  in  the  off-body  region  is  appropriate  provided  that  there  is  sufficient  grid  resolution  to  resolve  at  least  some 
of  the  turbulent  structures.  However,  despite  the  low  values  of  eddy  viscosity  in  the  wake  region  obtained  using  SA- 
DES  there  are  still  some  noticeable  qualitative  difference  between  using  SA-DES  and  inviscid  off-body  solvers.  In 
particular  the  wake  computed  using  the  SA-DES  turbulence  model  shows  fewer  small  scale  turbulent  structures  as 
well  as  a  narrower  overall  wake. 

The  convergence  histories  of  the  near-  and  off-body  solvers  are  plotted  in  figures  Figure  12(a)  and  Figure  12(b) 
respectively.  These  figures  show  that  solution  transients  are  decayed  after  25,000  time-steps  for  both  the  near-body 
and  off-body  solvers.  Figure  12(c)  depicts  the  computed  drag  coefficient  over  the  last  25,000  time-steps  of  the  time 
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(a)  Near-body  solver  (NSU3D)  convergence  history 


(b)  Off-body  solver  (SAMARC)  convergence  history 


(c)  Near-body  solver  (NSU3D)  Co  convergence 

Figure  12.  Vorticity  for  the  flow  over  a  sphere  with  a  SA-DES  turbulence  model  in  the  off-body  region  and  a  DDES  turbulence  model  in  the  near-body  region. 

history,  demonstrating  that  the  last  25,000  time-steps  are  suitable  for  obtaining  an  average  drag  coefficient  value. 
Therefore,  the  following  average  values  of  drag  are  computed  by  sampling  the  last  25,000  time-steps  of  the  solver. 

Table  1.  Computed  drag  coefficients  for  the  flow  over  a  sphere  with  flow  conditions  =  .3,  a  =  0°,  and  Re  =  6,760,000 


Euler-OB 

SA-RANS-OB 

SA-DES-OB 

Cd 

.411 

.550 

.457 

Table  1  contains  the  average  computed  drag  coefficients  for  the  flow  over  a  sphere.  The  computed  drag  coefficient 
obtained  utilizing  SA-RANS  in  the  off-body  mesh  is  significantly  higher  compared  with  those  computed  using  either 
of  the  inviscid  or  the  SA-DES  solvers  in  the  off-body  mesh.  Additionally,  it  is  interesting  to  note  that  using  SA-DES 
in  the  off-body  mesh  causes  the  drag  coefficient  to  increase  by  11.2%  over  the  drag  coefficient  computed  using  the 
inviscid  solver  in  the  off-body  mesh. 

This  test  case  has  demonstrated  that  the  effects  of  including  a  turbulence  model  in  the  off-body  mesh  can  be 
significant  to  both  qualitative  flow  features  as  well  as  relevant  simulation  outputs  such  as  drag  coefficient. 
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V.A.2.  Tilt  Rotor  Aeroacoustics  Model  (TRAM)  Rotor 

Helios  is  primarily  used  for  rotorcraft  aeromechanics  simulations.  In  order  to  examine  the  effects  of  including  turbu¬ 
lence  modeling  in  S  AMARC  on  a  revenant  rotorcraft  aeromechanics  problem,  this  work  considers  flow  though  the  Tilt 
Rotor  Aeroacoustics  Model  (TRAM)  rotor.  This  case  generates  strong  tip- vortices  that  must  be  adequately  resolved 
and  not  artificially  dissipated  by  the  off-body  solver.  Similarly  to  the  flow  over  a  sphere,  three  models  for  turbulence  in 
the  off-body  CFD  solver  are  compared:  inviscid,  SA-RANS,  and  SA-DES.  Qualitative  analysis  of  the  wake  structures 
is  discussed  to  examine  where  in  the  domain  and  how  the  SA-RANS  and  SA-DES  turbulence  models  are  altering  flow 
physics  .  Additionally,  effects  of  various  off-body  turbulence  models  on  the  figure  of  merit  (FM)  are  examined  as  a 
quantitative  basis  of  comparison. 

The  flow  conditions  of  the  TRAM  rotor  are  hover  conditions  with  a  tip  Mach  number  Mtip  =  0.625  and  collective 
pitch  0°  =  14°  on  the  rotor  blades.  The  tip  Reynolds  number  for  this  case  is  Retip  =  2.le6.  The  time-step  is  set  such 
that  the  rotor  rotates  .25°  degrees  per  time-step.  The  near-body  mesh  used  to  compute  this  flow  contains  9,273,348 
grid  points.  For  this  test  case  off-body  AMR  is  employed  every  10  time- steps  after  one  revolution  of  the  rotor  has  been 
completed.  For  all  test  cases  the  off-body  solver  is  initialized  with  a  uniform  cartesian  mesh  that  is  16  x  16  x  24.  The 
off-body  mesh  is  then  adapted  to  the  geometry  and  flow  physics  with  up  to  nine  levels  of  isotropic  mesh  refinement. 
Examples  of  the  types  of  meshes  used  in  these  computations  are  given  in  Figure  13(a)  and  Figure  13(b).  Using 


(a)  Computational  Mesh  at  x  =  0  plane  (b)  Computational  Mesh  at  x  =  0  plane 

Figure  13.  Flow  through  the  TRAM  rotor  in  hover,  using  the  SA  turbulence  model  in  the  off-body  mesh.  Vorticity  depicted  after  8  revolutions  of  the  rotor  and 
using  adaptive  mesh  refinement  in  the  off-body  region 


identical  inputs  three  variants  of  this  test  case  are  considered:  near-body  SA-RANS  with  off-body  inviscid,  near-body 
SA-RANS  with  off-body  SA-RANS,  and  near-body  SA-DDES  with  off-body  SA-DES. 

Figure  14(a)  depicts  the  vorticity  magnitude  contours  on  the  a  =  0  plane,  generated  by  using  an  inviscid  solver  in 
the  off-body  region.  This  figure  shows  that  the  tip  vortices  are  clearly  resolved  and  maintained  by  using  an  inviscid 
solver  in  the  off-body  region.  Furthermore,  there  is  little  evidence  of  strong  multi-scale  flow  structures  near  the  tip 
vortices.  However,  the  wake  generated  by  the  inboard  section  of  the  rotor  blades  as  well  as  by  the  hub  display  strong 
multi-scale  flow  physics  and  qualitatively  resemble  turbulence.  Figure  15(a)  depicts  the  vorticity  magnitude  contours 
on  the  a  =  0  plane,  generated  using  the  SA-RANS  turbulence  model  in  the  off-body  region.  One  should  immediate 
notice  from  Figure  15(a)  that  the  tip  vortices  are  also  clearly  resolved  despite  the  use  of  a  SA-RANS  turbulence 
model.  In  fact  comparing  Figure  14(a)  to  Figure  15(a)  shows  the  the  tip  vortex  structures  are  sufficiently  captured 
using  either  off-body  model.  Figure  15(b),  which  depicts  the  eddy  viscosity  for  this  case,  shows  why  the  tip  vortices 
are  not  dissipated  by  solver  as  the  eddy  viscosity  is  low  near  the  tip  vortices..  However,  examination  of  the  wake 
emanating  from  the  hub  shows  that  the  eddy  viscosity  if  relatively  large  in  this  region  and  therefore  the  wake  structure 
is  significantly  different  from  the  wake  structure  shown  in  Figure  14(a).  The  hub  wake  structure  in  Figure  15(a) 
contains  almost  no  small  scale  flow  structures  as  the  turbulence  model  is  modeling  all  length  scales  as  can  be  observed 
in  Figure  15(b).  Figure  16(a)  depicts  the  vorticity  magnitude  contours  on  the  a  =  0  plane,  generated  using  the  SA- 
DES  turbulence  model  in  the  off-body  region.  Similarly  to  the  previous  two  results,  Figure  15(a)  shows  that  the  tip 
vortices  are  clearly  resolved  and  comparable  to  the  inviscid  flow  solution  in  the  off-body  mesh.  The  SA-DES  model 
is  not  dissipating  the  vortices  because  as  with  the  SA-RANS  model,  the  eddy  viscosity  is  small  near  the  tip  vortices 
as  shown  in  Figure  16(b).  As  with  the  previous  two  off-body  flow  solutions  the  predominant  difference  between  the 
SA-DES,  SA-RANS,  and  inviscid  off-body  solutions  is  the  structure  of  the  wake  emanating  from  the  hub.  Examining 
Figure  16(a)  and  Figure  16(b)  shows  that  the  modeling  of  turbulence  is  relatively  active  in  the  hub  wake  when  using 
SA-DES  in  the  off-body  solution.  Furthermore,  comparing  Figure  14(a),  Figure  15(a)  and  Figure  16(a)  shows  that 
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Vorticity  Magnitude  [PLOT3D] 
o 


(a)  Vorticity  Contours  at  x  =  0  plane 


Figure  14.  Vorticity  magnitude  contours  for  flow  through  the  TRAM  rotor  in  hover,  with  an  inviscid  solver  in  the  off-body  mesh  after  8  revolutions  of  the  rotor 
and  using  adaptive  mesh  refinement  in  the  off-body  region 


(a)  Vorticity  Contours  at  x  =  0  plane  (b)  /jt  at  x  =  0  plane 


Figure  15.  Voritcity  and  eddy  viscosity  contours  on  the  x  =  0  plane  for  flow  through  the  TRAM  rotor  in  hover,  using  the  SA  turbulence  model  in  the  off-body 
mesh.  Solution  is  depicted  after  8  revolutions  of  the  rotor  and  using  adaptive  mesh  refinement  in  the  off-body  region 


Vorticity  Magnitude  [PLOT3D] 
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(a)  Vorticity  Contours  at  x  =  0  plane 

Figure  16.  Flow  through  the  TRAM  rotor  in  hover,  using  the  SA-DES  turbulence  model  in  the  off-body  mesh.  Vorticity  depicted  after  8  revolutions  of  the 
rotor  and  using  adaptive  mesh  refinement  in  the  off-body  region 


using  SA-DES  in  the  off-body  solver  results  in  a  hub  wake  that  is  somewhere  between  the  fully  modeled  SA-RANS 
solution  and  the  captured  only  inviscid  solution,  which  is  the  intended  and  expected  result. 
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(a)  FM  convergence  history  for  inviscid  off-body  solution  (b)  FM  convergence  history  for  SA-RANS  off-body  solution 


(c)  FM  convergence  history  for  RANS-DES  off-body  solution 

Figure  17.  Figure  of  merit  (FM)  convergence  history  for  the  flow  through  the  TRAM  rotor  using  inviscid,  SA-RANS,  and  SA-DES  in  the  off-body  region. 

Finally,  the  computed  figure  of  merit  (FM)  obtained  using  all  three  off-body  solution  options  is  compared.  Figure 
17(a)  through  Figure  17(c)  shows  the  convergence  histories  of  the  M’s  and  indicate  that  the  figure  of  merit  has 
become  converged  after  8  revolutions  of  the  rotor  are  completed.  Table  2  summarizes  the  average  FM  computed  using 
the  three  off-body  solutions  as  well  as  the  experimental  value.  Table  2  shows  that  at  the  spatial  and  temporal  resolution 


Table  2.  Computed  and  expire  mental  F M’s  for  the  flow  through  the  TRAM  rotor. 


Invsicd 

SA-RANS-OB 

SA-DES-OB 

Expriement 

FM 

.778 

.773 

.774 

.779 

considered  in  this  work  the  computed  figures  of  merit  are  approximately  equivalent  as  the  differences  between  all  three 
off-body  solver  options  are  very  small  (0(. 5%)).  Therefore,  this  work  has  validated  the  assumption  that  an  inviscid 
off-body  solver  is  an  appropriate  option  for  engineering  type  calculations. 

While  the  addition  of  turbulence  modeling  to  the  SAMARC  off-body  solver  has  not  improved  the  computed  FM, 
it  is  still  a  valuable  addition  SAMARC.  In  particular  has  some  sort  of  sub  grid  scale  model  is  beneficial  for  examine 
various  flow  features  that  can  occur  in  complex  rotorcraft  configurations. 

V.B.  Strand  Turbulence  Model  Validation 

The  S-A  implementation  in  the  strand  solver  is  validated  for  flow  over  a  flat  plate  at  M  =  0.2  and  Re  =  5  x  106,  based 
on  a  plate  of  length  unity.  This  case  was  taken  from  the  NASA  Langley  turbulence  modeling  resource  has  been  made 
for  this  case.23 
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The  grid  used  for  the  flat  plate  validation  is  a  137  x  97  grid  shown  in  Figure  18(a).  The  plate  leading  edge  begins 
at  x  =  0  and  extends  for  a  length  of  two.  A  short  inviscid  wall  entry  way  beginning  at  x  =  —0.33  is  provided  to  allow 
for  proper  inflow  conditions.  Stagnation  temperature  and  pressure  are  specified  at  the  inflow,  and  static  pressure  is 
specified  at  the  outflow.  The  turbulent  viscosity  field  for  this  case  is  shown  in  Figure  18(b),  which  has  been  scaled 
by  a  factor  of  40  vertically  to  facilitate  visualization.  Streamwise  velocity  and  turbulent  viscosity  profiles  are  shown 
in  Figure  19(a)  and  19(b)  for  two  locations  downstream  on  the  plate,  and  are  over  plotted  with  FUN3D  and  CFL3D 
results.  Note  that  good  agreement  is  obtained,  even  for  this  137  x  97  grid  which  is  16  times  more  coarse  than  the 
FUN3D  and  CFL3D  results  in  the  figures.  The  computed  drag  coefficient,  which  is  entirely  due  to  skin  friction  for  this 
case,  is  shown  in  Table  3,  along  with  FUN3D  and  CFL3D  results  for  the  same  grid.  The  drag  coefficient  falls  within 
the  range  predicted  by  the  established  codes. 

The  next  validation  case  is  2D  flow  over  a  NACA  0012  airfoil  at  M  =  0.15  and  Re  =  6  x  106  at  various  angles  of 
attack.  The  strand  grid  shown  in  Figure  20  consists  of  a  smoothed  NACA  0012  airfoil  containing  320  surface  nodes 
and  64  cells  along  each  strand,  yielding  20,480  total  strand  cells.  Figures  21  and  22  show  lift  coefficient  versus  angle 
of  attack  and  drag  coefficient  versus  lift  coefficient,  respectively,  along  with  the  corresponding  experimental  data  of 
Ladson.24  The  lift  data  is  matched  reasonably  well  for  all  angles  of  attack,  although  a  slight  over-prediction  of  lift 
is  observed  for  the  high-a  a  =  15°  case.  Drag  results  are  shown  in  Figure  22  with  the  a  =  15°  results  summarized 
in  Table  4.  FUN3D  and  CFL3D  results  using  a  much  finer  grid  (513  surface  nodes  instead  of  320)  are  included  for 
comparison.  In  general  the  strand  solver  results  fall  well  within  the  range  of  these  well  established  codes. 

VI.  Conclusions  &  Final  Paper 

This  paper  will  present  a  summary  of  turbulent  flow  validation  in  Helios  v3.  It  also  includes  a  validation  with  the 
new  strand  solver  that  is  eventually  intended  to  become  the  new  production  solver  in  Helios. 

Sample  calculations  are  shown  in  this  abstract,  the  final  paper  will  present  results  from  applications  as  defined  by 
the  AIAA  Fluid  Dynamics  Technical  Committee. 
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Figure  18.  Grid  and  turbulent  viscosity  contours  for  flow  over  a  flat  plate  at  M  =  0.2  and  Re  5  x  106. 
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(a)  x- velocity  profiles 


(b)  turbulent  viscosity  profile 


Figure  19.  Comparison  of  streamwise  velocity  and  turbulent  viscosity  profiles  for  flow  over  a  flat  plate  at  M  =  0.2  and  Re  5  x  106. 


Figure  20.  137  x  97  NACA  0012  strand  grid 
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Figure  21.  Q  vs.  a  compared  to  experiment  for  flow  over  NACA  0012  airfoil  at  M  =  0.15  and  Re  6  x  106. 
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Figure  22.  Q  vs.  Q  compared  to  experiment  for  flow  over  NACA  0012  airfoil  at  M  =  0.15  and  Re  6  x  106. 
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FUN3D  (quads) 

2.84005E-3 

FUN3D  (triangles) 

2.80289E-3 

CFL3D 

2.86621E-3 

Table  3.  Comparison  of  computed  drag  coefficients  for  flow  over  a  flat  plate  at  M  =  0.2  and  Re  =  5  x  106. 
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Cl 

cd 

Smoothed  strand  (320  points) 

1.5578 

0.02189 

FUN3D  (513  points) 

1.5547 

0.02159 

CFL3D  (513  points) 

1.5461 

0.02124 

Table  4.  Comparison  of  computed  lift  and  drag  coefficients  for  flow  over  a  NACA  0012  airfoil  at  M  —  0.15,  a  =  15°  Re  =  6  x  106.  Note  that  the  FUN3D  and 
CFL3D  results  use  a  much  finer  grid. 
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